R Code for GCB-19-0810
rm(list = ls()) # Remove all objects from memory
library(ggplot2)
library(gridExtra)
ENSO_colours<-c("blue", "red") # Type of anomaly
my_colours<-c("#feb24c93", "#fb6a4a93", "#9e9ac893", "#a6d96a93") # Thermal regime
theme_set (theme_classic() + theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
legend.position="none",
axis.text.x = element_text(angle = 90, vjust = 0.5),
plot.title = element_text(size=12, face="bold"),
# panel.border = element_blank(),
panel.border = element_rect(colour = "black", fill=NA, size=1)))
# Plot ENOS 1971 - 2014
ENOS<-read.table(file="http://www.cpc.ncep.noaa.gov/data/indices/ersst5.nino.mth.81-10.ascii",header=TRUE)
ENOS1982_2014 <-subset(ENOS, YR>1968 & YR<2015) # subset from 1970-2014
ENOS1982_2014$Date<-paste(ENOS1982_2014$YR, ENOS1982_2014$MON, sep = "-" )
ENOS1982_2014$Date<-paste(ENOS1982_2014$Date, "15", sep = "-" )
ENOS1982_2014$Date<-as.Date(ENOS1982_2014$Date, format = "%Y-%m-%d")
ENOS1982_2014$Anomaly <- factor(ifelse(ENOS1982_2014$ANOM.3 >=0, "Positive", "Negative"))
interp <- approx(ENOS1982_2014$Date, ENOS1982_2014$ANOM.3, n=10000)
DataENSO <- data.frame(Date=interp$x, Anomaly3=interp$y)
DataENSO$Anomaly <- factor(ifelse(DataENSO$Anomaly3>=0, "Positive", "Negative"))
ENSO<-ggplot(data=DataENSO, aes(x=Date, y=Anomaly3))+
scale_x_continuous("",expand=c(0, 0))+
scale_y_continuous("SST anomaly (ºC) in Niño.3",expand=c(0, 0))+
geom_area(aes(fill=Anomaly), alpha=0.7)+
scale_fill_manual(values=ENSO_colours) + labs(title="a")+
geom_hline(yintercept=0, linetype="solid", color="black")+
geom_line(position = "identity", size=0.3)+
theme(axis.text.x=element_blank(),
axis.title.x=element_blank(),
axis.ticks.x = element_blank())+
annotate("text", x = 1900, y = 2.0, label = "El Niño \n72-73", size=3)+
annotate("text", x = 5600, y = 2.0, label = "El Niño \n82-83", size=3)+
annotate("text", x = 11100, y = 2.0, label = "El Niño \n97-98", size=3)
ENSO
# ENSOb<-ggplot(data=ENOS1982_2014, aes(x=Date, y=ANOM.3))+
# geom_hline(yintercept=0, linetype="solid", color="black")+
# scale_x_date("", limit=c(as.Date("1969-01-15"), as.Date("2014-12-31")),
# breaks= "12 months",
# expand=c(0, 0))+
# scale_y_continuous("SST anomaly (ºC) in Niño.3",expand=c(0, 0))+
# geom_area(aes(fill=Anomaly))+
# geom_line(position = "identity")+
# scale_fill_manual(values=ENSO_colours) + labs(title="(a)")
# ENSOb
# Select data
cover <- read.csv("Data/Coral_Cover_ETP.csv", header=TRUE) # Select file: Coral_Cover_ETP.csv
# Convert to numeric variable the response variable and convert to nominal variables the factors
cover$Coral_cover <- as.numeric(cover$Coral_cover) # convert to numeric variable
cover$Region <- as.factor(cover$Region) # convert to nominal factors
cover$Location <- as.factor(cover$Location)
cover$Site <- as.factor(cover$Site)
cover$Country <- as.factor(cover$Country)
cover$Thermal_regime <- factor(cover$Thermal_regime,
levels=c("Thermally_Stable", "Upwelling", "Galapagos_Islands","Seasonal"),
labels=c("Thermally stable", "Tropical upwelling", "Equatorial upwelling","Seasonal"))
cover$Lat <- as.numeric(cover$Lat)
cover$Lon <- as.numeric(cover$Lon)
cover$Year <- as.numeric(cover$Year)
#cover$Year <- as.numeric((cover$Year)-1970) # transform years to begin in 0
summary(cover)
## ID1 Reference Province
## Min. : 1.0 Fong, 2017 : 33 Cortezian : 10
## 1st Qu.:144.2 Glynn, 1983 : 31 Galapagos : 61
## Median :285.5 Calderón-Aguilera, 2006: 26 Tropical_East_Pacific:495
## Mean :285.1 Glynn, 1990 : 24
## 3rd Qu.:426.8 Guzman, 2004 : 24
## Max. :568.0 Cortés, 2010 : 23
## (Other) :405
## Ecoregion Region
## Nicoya :238 Gulf_of_Chiriqui :133
## Panama_Bight :117 Colombia : 75
## Cortezian : 71 Mexican_Central : 72
## Mexican_Tropical_Pacific : 50 Costa_Rica_Southern : 56
## Eastern_Galapagos_Islands : 31 Nicaragua_Papagayo_Zone: 49
## Northern_Galapagos_Islands: 28 Panama_Bight : 42
## (Other) : 31 (Other) :139
## Country Platform Thermal_regime
## Panama :151 Continental:470 Thermally stable :250
## Costa_Rica:143 Oceanic : 96 Tropical upwelling :177
## Mexico :123 Equatorial upwelling: 61
## Colombia : 75 Seasonal : 78
## Ecuador : 61
## Nicaragua : 11
## (Other) : 2
## Location Site Year Date
## Gulf_of_Chiriqui:121 Uva_Island : 67 Min. :1970 1/1/09 : 55
## Gorgona_Island : 56 La_Azufrada : 29 1st Qu.:1987 1/1/06 : 32
## Pearl_Islands : 40 Cano_island : 27 Median :1999 1/1/02 : 31
## Bahia_Banderas : 33 Gorgona_Island: 26 Mean :1997 1/1/03 : 31
## Gulf_of_Papagayo: 30 Saboga_Island : 24 3rd Qu.:2006 1/1/97 : 26
## Los_Cabos : 28 Cabo_Pulmo : 17 Max. :2014 1/1/98 : 22
## (Other) :258 (Other) :376 (Other):369
## Period_ENOS Coral_cover Algal_cover Lat
## 1971-1972: 24 Min. : 0.000 Min. : 0.00 Min. :-1.352
## 1973-1982: 60 1st Qu.: 8.568 1st Qu.:12.61 1st Qu.: 5.551
## 1983-1997:169 Median :25.000 Median :35.84 Median : 8.578
## 1998-2016:313 Mean :29.185 Mean :37.57 Mean : 9.499
## 3rd Qu.:45.725 3rd Qu.:54.81 3rd Qu.:10.791
## Max. :95.720 Max. :98.00 Max. :25.796
## NA's :390
## Lon
## Min. :-111.23
## 1st Qu.: -91.82
## Median : -83.87
## Mean : -87.92
## 3rd Qu.: -81.75
## Max. : -77.34
##
# Aggregate by site
aggr.site <- aggregate(Coral_cover ~ Year+Region+Location+Site+Period_ENOS+Thermal_regime,FUN=mean,data=cover)
# Aggregate by location
aggr.location <- aggregate(Coral_cover~Year+Region+Location+Period_ENOS+Thermal_regime,FUN=mean,data=aggr.site)
# Aggregate by region
aggr.region <- aggregate(Coral_cover~Year+Region+Period_ENOS+Thermal_regime,FUN=mean,data=aggr.location)
# Create groups for each of three ENSO time-intervals
ENSO.73_82 <- subset(aggr.region, Period_ENOS == "1973-1982")
ENSO.83_97 <- subset(aggr.region, Period_ENOS == "1983-1997")
ENSO.98_14 <- subset(aggr.region, Period_ENOS == "1998-2016")
# Exclude regions with cero and only one observation of coral cover per time-interval
ENSO.83_97<- ENSO.83_97[!ENSO.83_97$Region == "Costa_Rica_Central", ]
ENSO.83_97<- ENSO.83_97[!ENSO.83_97$Region == "El_Salvador", ]
ENSO.83_97<- ENSO.83_97[!ENSO.83_97$Region == "Western_Galapagos_Islands", ]
ENSO.83_97<- ENSO.83_97[!ENSO.83_97$Region == "Clipperton", ]
ENSO.83_97<- ENSO.83_97[!ENSO.83_97$Region == "Mexican_Southern", ]
ENSO.83_97<- ENSO.83_97[!ENSO.83_97$Region == "Northern_Galapagos_Islands", ]
ENSO.98_14<- ENSO.98_14[!ENSO.98_14$Region == "Clipperton", ]
ENSO.98_14<- ENSO.98_14[!ENSO.98_14$Region == "El_Salvador", ]
ENSO.98_14<- ENSO.98_14[!ENSO.98_14$Region == "Western_Galapagos_Islands", ]
ENSO.98_14<- ENSO.98_14[!ENSO.98_14$Region == "Revillagigedo_Islands", ]
ENSO.83_14 <- rbind(ENSO.83_97,ENSO.98_14) # Second and third time intervals
ENSO.73_14 <- rbind(ENSO.73_82,ENSO.83_14) # First, second and third time intervals
ENSO.73_97 <- rbind(ENSO.73_82,ENSO.83_97) # First and second and time intervals
library(nlme)
# Coral_cover = response variable
# Region = repeated measure "subject" or repeated experimental unit (random factor)
# Year = continuous numeric variable
# Thermal_regime = fixed factor
# All years 1970_2014
model1 <- lme(
Coral_cover ~ -1 + Thermal_regime * Year, random = ~1|Region, data=aggr.region)
summary(model1)
## Linear mixed-effects model fit by REML
## Data: aggr.region
## AIC BIC logLik
## 1727.902 1760.58 -853.9508
##
## Random effects:
## Formula: ~1 | Region
## (Intercept) Residual
## StdDev: 12.5915 16.52254
##
## Fixed effects: Coral_cover ~ -1 + Thermal_regime * Year
## Value Std.Error DF t-value
## Thermal_regimeThermally stable -428.2896 377.6124 181 -1.1342042
## Thermal_regimeTropical upwelling -323.5241 371.8720 14 -0.8699878
## Thermal_regimeEquatorial upwelling 254.5442 573.7481 14 0.4436514
## Thermal_regimeSeasonal 1036.7134 843.3488 181 1.2292819
## Year 0.2264 0.1889 181 1.1987453
## Thermal_regimeTropical upwelling:Year -0.0478 0.2650 181 -0.1804381
## Thermal_regimeEquatorial upwelling:Year -0.3448 0.3436 181 -1.0033294
## Thermal_regimeSeasonal:Year -0.7293 0.4592 181 -1.5880369
## p-value
## Thermal_regimeThermally stable 0.2582
## Thermal_regimeTropical upwelling 0.3990
## Thermal_regimeEquatorial upwelling 0.6641
## Thermal_regimeSeasonal 0.2206
## Year 0.2322
## Thermal_regimeTropical upwelling:Year 0.8570
## Thermal_regimeEquatorial upwelling:Year 0.3170
## Thermal_regimeSeasonal:Year 0.1140
## Correlation:
## Thr_Ts Thr_Tu Thr_Eu Thrm_S Year
## Thermal_regimeTropical upwelling 0.000
## Thermal_regimeEquatorial upwelling 0.000 0.000
## Thermal_regimeSeasonal 0.014 0.000 0.000
## Year -1.000 0.000 0.000 -0.015
## Thermal_regimeTropical upwelling:Year 0.713 -0.701 0.000 0.010 -0.713
## Thermal_regimeEquatorial upwelling:Year 0.550 0.000 -0.835 0.008 -0.550
## Thermal_regimeSeasonal:Year 0.398 0.000 0.000 -0.912 -0.397
## T_Tu:Y T_Eu:Y
## Thermal_regimeTropical upwelling
## Thermal_regimeEquatorial upwelling
## Thermal_regimeSeasonal
## Year
## Thermal_regimeTropical upwelling:Year
## Thermal_regimeEquatorial upwelling:Year 0.392
## Thermal_regimeSeasonal:Year 0.283 0.218
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.323212218 -0.685668507 -0.007907488 0.578042210 3.071703318
##
## Number of Observations: 202
## Number of Groups: 16
anova(model1)
## numDF denDF F-value p-value
## Thermal_regime 4 179 14.698982 <.0001
## Year 1 179 0.717103 0.3982
## Thermal_regime:Year 3 179 1.093597 0.3532
plot(ranef(model1)) # shows symmetrical scatter effects around zero
plot(model1) # plot residuals vs fitted
resnorm1 <- resid(model1)
hist(resnorm1, xlab = "Residuals", main = "") # residuals are normally distributed?
coef.m1 <- as.data.frame(coef(summary(model1))) # Coefficients of the model
plot(model1, resid(., type = "p") ~ fitted(.) | Region, abline = 0)
plot(model1, Coral_cover ~ fitted(.) | Region, abline = c(0,1))
plot(model1, Coral_cover ~ fitted(.) | Thermal_regime, abline = c(0,1))
plot(model1, Coral_cover ~ fitted(.), abline = c(0,1))
#Using raw data (Only thermal regime is significant in model)
Model1_plot <- ggplot(aggr.region, aes(x=Year, y=Coral_cover)) +
geom_boxplot(aes(group=Year), outlier.shape = NA) +
stat_boxplot(aes(group=Year), geom = 'errorbar')+
geom_smooth(method = lm, se=FALSE, linetype = "dashed", colour="deeppink")+
geom_smooth(span = 0.2, se=T, colour="aquamarine")+
#geom_point(aes(colour=Thermal_regime), alpha=0.5) + ylim(0,79) +
scale_y_continuous("Live coral cover (%)", limits = c(-2, 80), expand = c(0,0))+
scale_x_continuous("", limits = c(1969, 2015),
breaks = seq(1970, 2014, by=5), expand = c(0,0))+
scale_color_manual(values=my_colours) + labs(title="b")+
#annotate("text", x = 1983, y = 70, label = "*", size=8)+
annotate("rect", xmin = 1982, xmax = 1983, ymin = 0, ymax = 80, alpha = .2, fill="red")+
annotate("rect", xmin = 1997, xmax = 1998, ymin = 0, ymax = 80, alpha = .2, fill="red")+
annotate("rect", xmin = 1972, xmax = 1973, ymin = 0, ymax = 80, alpha = .2, fill="red")
Model1_plot
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
# Time-interval 1 - ENSO.73_82
model2 <- lme(Coral_cover ~ -1 + Thermal_regime * Year, random = ~1|Region, data=ENSO.73_82)
summary(model2)
## Linear mixed-effects model fit by REML
## Data: ENSO.73_82
## AIC BIC logLik
## 164.3022 170.4829 -74.15108
##
## Random effects:
## Formula: ~1 | Region
## (Intercept) Residual
## StdDev: 20.41219 11.23257
##
## Fixed effects: Coral_cover ~ -1 + Thermal_regime * Year
## Value Std.Error DF t-value
## Thermal_regimeThermally stable -570.074 2795.777 4 -0.2039053
## Thermal_regimeTropical upwelling -8303.810 3542.253 4 -2.3442170
## Thermal_regimeEquatorial upwelling -3307.477 3419.771 4 -0.9671634
## Year 0.301 1.412 13 0.2131857
## Thermal_regimeTropical upwelling:Year 3.912 2.281 13 1.7150859
## Thermal_regimeEquatorial upwelling:Year 1.381 2.232 13 0.6189979
## p-value
## Thermal_regimeThermally stable 0.8484
## Thermal_regimeTropical upwelling 0.0790
## Thermal_regimeEquatorial upwelling 0.3882
## Year 0.8345
## Thermal_regimeTropical upwelling:Year 0.1101
## Thermal_regimeEquatorial upwelling:Year 0.5466
## Correlation:
## Thr_Ts Thr_Tu Thr_Eu Year T_Tu:Y
## Thermal_regimeTropical upwelling 0.000
## Thermal_regimeEquatorial upwelling 0.000 0.000
## Year -1.000 0.000 0.000
## Thermal_regimeTropical upwelling:Year 0.619 -0.785 0.000 -0.619
## Thermal_regimeEquatorial upwelling:Year 0.633 0.000 -0.774 -0.633 0.392
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.4650132 -0.4802879 -0.2103293 0.5465022 1.8906401
##
## Number of Observations: 22
## Number of Groups: 7
anova(model2)
## numDF denDF F-value p-value
## Thermal_regime 3 4 3.434033 0.1322
## Year 1 13 3.583180 0.0808
## Thermal_regime:Year 2 13 1.472433 0.2652
plot(ranef(model2))
plot(model2)
resnorm2 <- resid(model2)
hist(resnorm2, xlab = "Residuals", main = "")
coef.m2 <- as.data.frame(coef(summary(model2)))
plot(model2, resid(., type = "p") ~ fitted(.) | Region, abline = 0)
plot(model2, Coral_cover ~ fitted(.) | Region, abline = c(0,1))
plot(model2, Coral_cover ~ fitted(.) | Thermal_regime, abline = c(0,1))
# Note that this first time interval has low number of regions at each year, even the
# GLMM model is adjusted, please warning with interpretation/generalization of results
#Using raw data (the model is not significant)
Model2_plot <- ggplot(ENSO.73_82, aes(x=Year, y=Coral_cover, colour=Thermal_regime)) +
geom_smooth(method = lm, se=FALSE, linetype = "dashed", alpha=1, size=2)+
geom_point(aes(fill=factor(Thermal_regime)),
shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.5) +
scale_y_continuous("Live coral cover (%)", limits = c(0, 80))+
scale_x_continuous("", limits = c(1973.7, 1982.3),
breaks = seq(1974, 1982, by=2), expand = c(0.02,0.02))+
scale_fill_manual(values=my_colours) +
scale_colour_manual(values=my_colours)+
labs(title="c")
Model2_plot
# Time-interval 3 - ENSO.83_97
model3 <- lme(Coral_cover ~ -1 + Thermal_regime * Year, random = ~1|Region, data=ENSO.83_97)
summary(model3)
## Linear mixed-effects model fit by REML
## Data: ENSO.83_97
## AIC BIC logLik
## 476.7179 497.1484 -228.3589
##
## Random effects:
## Formula: ~1 | Region
## (Intercept) Residual
## StdDev: 4.741945 9.686494
##
## Fixed effects: Coral_cover ~ -1 + Thermal_regime * Year
## Value Std.Error DF t-value
## Thermal_regimeThermally stable -2536.003 926.2417 50 -2.737950
## Thermal_regimeTropical upwelling -7100.953 882.5060 8 -8.046351
## Thermal_regimeEquatorial upwelling 111.787 2814.6054 8 0.039717
## Thermal_regimeSeasonal 6083.292 1940.6223 50 3.134712
## Year 1.281 0.4654 50 2.753386
## Thermal_regimeTropical upwelling:Year 2.306 0.6426 50 3.589323
## Thermal_regimeEquatorial upwelling:Year -1.337 1.4905 50 -0.897275
## Thermal_regimeSeasonal:Year -4.319 1.0734 50 -4.023670
## p-value
## Thermal_regimeThermally stable 0.0085
## Thermal_regimeTropical upwelling 0.0000
## Thermal_regimeEquatorial upwelling 0.9693
## Thermal_regimeSeasonal 0.0029
## Year 0.0082
## Thermal_regimeTropical upwelling:Year 0.0008
## Thermal_regimeEquatorial upwelling:Year 0.3739
## Thermal_regimeSeasonal:Year 0.0002
## Correlation:
## Thr_Ts Thr_Tu Thr_Eu Thrm_S Year
## Thermal_regimeTropical upwelling 0.000
## Thermal_regimeEquatorial upwelling 0.000 0.000
## Thermal_regimeSeasonal 0.014 0.000 0.000
## Year -1.000 0.000 0.000 -0.014
## Thermal_regimeTropical upwelling:Year 0.724 -0.689 0.000 0.010 -0.724
## Thermal_regimeEquatorial upwelling:Year 0.312 0.000 -0.950 0.005 -0.312
## Thermal_regimeSeasonal:Year 0.420 0.000 0.000 -0.901 -0.420
## T_Tu:Y T_Eu:Y
## Thermal_regimeTropical upwelling
## Thermal_regimeEquatorial upwelling
## Thermal_regimeSeasonal
## Year
## Thermal_regimeTropical upwelling:Year
## Thermal_regimeEquatorial upwelling:Year 0.226
## Thermal_regimeSeasonal:Year 0.304 0.131
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.99663110 -0.53542294 -0.06265663 0.58765352 2.99408354
##
## Number of Observations: 65
## Number of Groups: 10
anova(model3)
## numDF denDF F-value p-value
## Thermal_regime 4 48 62.03136 <.0001
## Year 1 48 39.45390 <.0001
## Thermal_regime:Year 3 48 14.56275 <.0001
plot(ranef(model3))
plot(model3)
resnorm3 <- resid(model3)
hist(resnorm3, xlab = "Residuals", main = "")
coef.m3 <- as.data.frame(coef(summary(model3)))
plot(model3, resid(., type = "p") ~ fitted(.) | Region, abline = 0)
plot(model3, Coral_cover ~ fitted(.) | Region, abline = c(0,1))
plot(model3, Coral_cover ~ fitted(.) | Thermal_regime, abline = c(0,1))
# Create a new data frame for independent variables
NewData_3 <- expand.grid(Thermal_regime=unique(ENSO.83_97$Thermal_regime),
#Region=unique(aggr.region$Region),
Year=seq((min(ENSO.83_97$Year)), (max(ENSO.83_97$Year))))
pred_3 <- predict(model3, newdata=NewData_3, level=0, asList = FALSE)
# Using model data
Model3b_plot <- ggplot(ENSO.83_97, aes(x=Year, y=Coral_cover, colour=Thermal_regime)) +
geom_point(aes(fill=factor(Thermal_regime)),
shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.5) +
scale_y_continuous(" ", limits = c(0, 80))+
scale_x_continuous("", limits = c(1982.7, 1997.3),
breaks = seq(1983, 1997, by=2), expand = c(0.02,0.02))+
geom_line(data=NewData_3,
aes(y=predict(model3, level=0, newdata=NewData_3)), size=2) +
theme(axis.text.y=element_blank(),
axis.title.y=element_blank(),
legend.position = c(0.58, -0.4),
legend.title=element_blank(),
legend.text=element_text(size=6),
legend.box ="horizontal")+
guides(col = guide_legend(nrow = 2))+
scale_fill_manual(values=my_colours) +
scale_colour_manual(values=my_colours)+
labs(title="d")
Model3b_plot
# Time-interval 4 - ENSO.98_14
model4 <- lme(Coral_cover ~ -1 + Thermal_regime * Year, random = ~1|Region, data=ENSO.98_14)
summary(model4)
## Linear mixed-effects model fit by REML
## Data: ENSO.98_14
## AIC BIC logLik
## 866.8707 892.5142 -423.4354
##
## Random effects:
## Formula: ~1 | Region
## (Intercept) Residual
## StdDev: 11.83418 15.29875
##
## Fixed effects: Coral_cover ~ -1 + Thermal_regime * Year
## Value Std.Error DF t-value
## Thermal_regimeThermally stable -1881.9759 1396.7590 87 -1.3473877
## Thermal_regimeTropical upwelling 2994.9557 1056.7535 10 2.8341102
## Thermal_regimeEquatorial upwelling -1881.0600 2875.6287 10 -0.6541387
## Thermal_regimeSeasonal 247.4085 1577.1131 87 0.1568743
## Year 0.9516 0.6965 87 1.3661224
## Thermal_regimeTropical upwelling:Year -2.4271 0.8732 87 -2.7795278
## Thermal_regimeEquatorial upwelling:Year -0.0018 1.5937 87 -0.0011048
## Thermal_regimeSeasonal:Year -1.0581 1.0458 87 -1.0117420
## p-value
## Thermal_regimeThermally stable 0.1814
## Thermal_regimeTropical upwelling 0.0177
## Thermal_regimeEquatorial upwelling 0.5278
## Thermal_regimeSeasonal 0.8757
## Year 0.1754
## Thermal_regimeTropical upwelling:Year 0.0067
## Thermal_regimeEquatorial upwelling:Year 0.9991
## Thermal_regimeSeasonal:Year 0.3145
## Correlation:
## Thr_Ts Thr_Tu Thr_Eu Thrm_S Year
## Thermal_regimeTropical upwelling 0.000
## Thermal_regimeEquatorial upwelling 0.000 0.000
## Thermal_regimeSeasonal 0.008 0.000 0.000
## Year -1.000 0.000 0.000 -0.008
## Thermal_regimeTropical upwelling:Year 0.798 -0.603 0.000 0.006 -0.798
## Thermal_regimeEquatorial upwelling:Year 0.437 0.000 -0.899 0.003 -0.437
## Thermal_regimeSeasonal:Year 0.660 0.000 0.000 -0.746 -0.660
## T_Tu:Y T_Eu:Y
## Thermal_regimeTropical upwelling
## Thermal_regimeEquatorial upwelling
## Thermal_regimeSeasonal
## Year
## Thermal_regimeTropical upwelling:Year
## Thermal_regimeEquatorial upwelling:Year 0.349
## Thermal_regimeSeasonal:Year 0.526 0.288
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.98327607 -0.50110676 -0.02557425 0.30976474 2.71954447
##
## Number of Observations: 104
## Number of Groups: 12
anova(model4)
## numDF denDF F-value p-value
## Thermal_regime 4 85 15.541089 <.0001
## Year 1 85 1.227329 0.2711
## Thermal_regime:Year 3 85 2.983540 0.0358
plot(ranef(model4))
plot(model4)
resnorm4 <- resid(model4)
hist(resnorm4, xlab = "Residuals", main = "")
coef.m4 <- as.data.frame(coef(summary(model4)))
plot(model4, resid(., type = "p") ~ fitted(.) | Region, abline = 0)
plot(model4, Coral_cover ~ fitted(.) | Region, abline = c(0,1))
plot(model4, Coral_cover ~ fitted(.) | Thermal_regime, abline = c(0,1))
# Competing model with autocorrelation structure
model4b <- lme(Coral_cover ~ -1 + Thermal_regime * Year, random = ~1|Region,
data=ENSO.98_14, correlation=corCompSymm(form=~Year|Region))
anova(model4, model4b) # ~ same results
## Model df AIC BIC logLik Test L.Ratio p-value
## model4 1 10 866.8707 892.5142 -423.4354
## model4b 2 11 868.8707 897.0786 -423.4354 1 vs 2 4.547474e-13 1
# Create a new data frame for independent variables
NewData_4 <- expand.grid(Thermal_regime=unique(ENSO.98_14$Thermal_regime),
#Region=unique(aggr.region$Region),
Year=seq((min(ENSO.98_14$Year)), (max(ENSO.98_14$Year))))
pred_4 <- predict(model4, newdata=NewData_4, level=0)
#summary(pred)
#length(pred)
# Using model data
Model4b_plot <- ggplot(ENSO.98_14, aes(x=Year, y=Coral_cover, colour=Thermal_regime)) +
geom_point(aes(fill=factor(Thermal_regime)),
shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.5) +
scale_y_continuous(" ", limits = c(0, 80))+
scale_x_continuous("", limits = c(1997.7, 2014.3),
breaks = seq(1998, 2014, by=2), expand = c(0.02,0.02))+
geom_line(data=NewData_4, aes(y=predict(model4, level=0, newdata=NewData_4)), size=2)+
theme(axis.text.y=element_blank(), axis.title.y=element_blank())+
scale_fill_manual(values=my_colours) +
scale_colour_manual(values=my_colours)+
labs(title="e")
Model4b_plot
# Time-interval 5 - ENSO.83_14
model5 <- lme(Coral_cover ~ -1 + Thermal_regime * Year, random = ~1|Region, data=ENSO.83_14)
summary(model5)
## Linear mixed-effects model fit by REML
## Data: ENSO.83_14
## AIC BIC logLik
## 1423.169 1453.983 -701.5844
##
## Random effects:
## Formula: ~1 | Region
## (Intercept) Residual
## StdDev: 12.53351 15.47206
##
## Fixed effects: Coral_cover ~ -1 + Thermal_regime * Year
## Value Std.Error DF t-value
## Thermal_regimeThermally stable -1736.8288 514.4829 151 -3.375873
## Thermal_regimeTropical upwelling 36.1048 507.2122 11 0.071183
## Thermal_regimeEquatorial upwelling -1334.0770 907.7813 11 -1.469602
## Thermal_regimeSeasonal 1031.1175 792.1371 151 1.301691
## Year 0.8796 0.2572 151 3.420336
## Thermal_regimeTropical upwelling:Year -0.8804 0.3609 151 -2.439274
## Thermal_regimeEquatorial upwelling:Year -0.2024 0.5212 151 -0.388367
## Thermal_regimeSeasonal:Year -1.3793 0.4685 151 -2.943936
## p-value
## Thermal_regimeThermally stable 0.0009
## Thermal_regimeTropical upwelling 0.9445
## Thermal_regimeEquatorial upwelling 0.1697
## Thermal_regimeSeasonal 0.1950
## Year 0.0008
## Thermal_regimeTropical upwelling:Year 0.0159
## Thermal_regimeEquatorial upwelling:Year 0.6983
## Thermal_regimeSeasonal:Year 0.0038
## Correlation:
## Thr_Ts Thr_Tu Thr_Eu Thrm_S Year
## Thermal_regimeTropical upwelling 0.000
## Thermal_regimeEquatorial upwelling 0.000 0.000
## Thermal_regimeSeasonal 0.015 0.000 0.000
## Year -1.000 0.000 0.000 -0.016
## Thermal_regimeTropical upwelling:Year 0.712 -0.702 0.000 0.011 -0.712
## Thermal_regimeEquatorial upwelling:Year 0.493 0.000 -0.870 0.008 -0.493
## Thermal_regimeSeasonal:Year 0.536 0.000 0.000 -0.836 -0.535
## T_Tu:Y T_Eu:Y
## Thermal_regimeTropical upwelling
## Thermal_regimeEquatorial upwelling
## Thermal_regimeSeasonal
## Year
## Thermal_regimeTropical upwelling:Year
## Thermal_regimeEquatorial upwelling:Year 0.352
## Thermal_regimeSeasonal:Year 0.381 0.264
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.52794002 -0.47964749 -0.06607999 0.50350339 3.04770181
##
## Number of Observations: 169
## Number of Groups: 13
anova(model5)
## numDF denDF F-value p-value
## Thermal_regime 4 149 14.198648 <.0001
## Year 1 149 4.270123 0.0405
## Thermal_regime:Year 3 149 3.799303 0.0116
plot(ranef(model5))
plot(model5)
resnorm5 <- resid(model5)
hist(resnorm5, xlab = "Residuals", main = "")
coef.m5 <- as.data.frame(coef(summary(model5)))
plot(model5, resid(., type = "p") ~ fitted(.) | Region, abline = 0)
plot(model5, Coral_cover ~ fitted(.) | Region, abline = c(0,1))
plot(model5, Coral_cover ~ fitted(.) | Thermal_regime, abline = c(0,1))
# Create a new data frame for independent variables
NewData_5 <- expand.grid(Thermal_regime=unique(ENSO.83_14$Thermal_regime),
#Region=unique(aggr.region$Region),
Year=seq((min(ENSO.83_14$Year)), (max(ENSO.83_14$Year))))
pred_5 <- predict(model5, newdata=NewData_5, level=0)
# Using model data
Model5b_plot <- ggplot(ENSO.83_14, aes(x=Year, y=Coral_cover, colour=Thermal_regime)) +
geom_point(aes(fill=factor(Thermal_regime)),
shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.5) +
scale_y_continuous("Live coral cover (%)", limits = c(0, 80))+
scale_x_continuous("", limits = c(1982.7, 2014.3),
breaks = seq(1983, 2014, by=2), expand = c(0.02,0.02))+
geom_line(data=NewData_5, aes(y=predict(model5, level=0, newdata=NewData_5)), size=2)+
theme(axis.text.y=element_blank(),
axis.title.y=element_blank(),
legend.position = c(0.58, -0.3),
legend.title=element_blank(),
legend.text=element_text(size=6),
legend.direction="horizontal")+
scale_fill_manual(values=my_colours) +
scale_colour_manual(values=my_colours)+
labs(title="Model 5")
Model5b_plot
library(multcomp)
ENSO.73_97$Year <- as.factor(ENSO.73_97$Year)
comp_years_ENSO.73_97 <- aov(Coral_cover ~ Year + Thermal_regime, data = ENSO.73_97)
m.comp_ENSO.73_97 <- glht(comp_years_ENSO.73_97, linfct = mcp(Year = "Tukey"))
summary(m.comp_ENSO.73_97, test = univariate())
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = Coral_cover ~ Year + Thermal_regime, data = ENSO.73_97)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 1975 - 1974 == 0 -3.44519 12.08474 -0.285 0.776577
## 1976 - 1974 == 0 11.67379 13.58569 0.859 0.393669
## 1977 - 1974 == 0 2.20768 17.29237 0.128 0.898846
## 1978 - 1974 == 0 -21.04773 13.64079 -1.543 0.128178
## 1979 - 1974 == 0 12.32727 13.64079 0.904 0.369826
## 1980 - 1974 == 0 2.61704 12.23990 0.214 0.831431
## 1981 - 1974 == 0 10.89268 17.29237 0.630 0.531184
## 1982 - 1974 == 0 9.29041 10.82386 0.858 0.394186
## 1983 - 1974 == 0 -26.56453 12.25535 -2.168 0.034237 *
## 1984 - 1974 == 0 -17.28488 11.32342 -1.526 0.132236
## 1985 - 1974 == 0 -6.66044 13.76573 -0.484 0.630289
## 1986 - 1974 == 0 -16.73332 11.32342 -1.478 0.144790
## 1987 - 1974 == 0 -6.97546 10.29300 -0.678 0.500618
## 1988 - 1974 == 0 -5.21764 11.53808 -0.452 0.652776
## 1989 - 1974 == 0 -9.53355 10.53984 -0.905 0.369396
## 1990 - 1974 == 0 -10.69232 17.29237 -0.618 0.538741
## 1991 - 1974 == 0 4.34819 12.74709 0.341 0.734231
## 1992 - 1974 == 0 0.52130 10.65659 0.049 0.961150
## 1993 - 1974 == 0 -6.10824 11.04278 -0.553 0.582255
## 1994 - 1974 == 0 3.00209 11.00761 0.273 0.786014
## 1995 - 1974 == 0 25.24915 11.51467 2.193 0.032279 *
## 1996 - 1974 == 0 8.32457 11.51467 0.723 0.472565
## 1997 - 1974 == 0 0.07463 10.75923 0.007 0.994489
## 1998 - 1974 == 0 26.29686 17.32517 1.518 0.134395
## 1976 - 1975 == 0 15.11898 13.58569 1.113 0.270282
## 1977 - 1975 == 0 5.65287 17.29237 0.327 0.744901
## 1978 - 1975 == 0 -17.60254 13.64079 -1.290 0.201934
## 1979 - 1975 == 0 15.77246 13.64079 1.156 0.252231
## 1980 - 1975 == 0 6.06223 12.23990 0.495 0.622239
## 1981 - 1975 == 0 14.33787 17.29237 0.829 0.410365
## 1982 - 1975 == 0 12.73559 10.82386 1.177 0.244072
## 1983 - 1975 == 0 -23.11934 12.25535 -1.886 0.064156 .
## 1984 - 1975 == 0 -13.83970 11.32342 -1.222 0.226485
## 1985 - 1975 == 0 -3.21526 13.76573 -0.234 0.816128
## 1986 - 1975 == 0 -13.28814 11.32342 -1.174 0.245307
## 1987 - 1975 == 0 -3.53027 10.29300 -0.343 0.732834
## 1988 - 1975 == 0 -1.77246 11.53808 -0.154 0.878435
## 1989 - 1975 == 0 -6.08837 10.53984 -0.578 0.565697
## 1990 - 1975 == 0 -7.24713 17.29237 -0.419 0.676669
## 1991 - 1975 == 0 7.79338 12.74709 0.611 0.543293
## 1992 - 1975 == 0 3.96649 10.65659 0.372 0.711070
## 1993 - 1975 == 0 -2.66305 11.04278 -0.241 0.810269
## 1994 - 1975 == 0 6.44728 11.00761 0.586 0.560304
## 1995 - 1975 == 0 28.69434 11.51467 2.492 0.015533 *
## 1996 - 1975 == 0 11.76975 11.51467 1.022 0.310882
## 1997 - 1975 == 0 3.51982 10.75923 0.327 0.744717
## 1998 - 1975 == 0 29.74205 17.32517 1.717 0.091281 .
## 1977 - 1976 == 0 -9.46611 18.35985 -0.516 0.608069
## 1978 - 1976 == 0 -32.72152 15.11308 -2.165 0.034435 *
## 1979 - 1976 == 0 0.65348 15.11308 0.043 0.965657
## 1980 - 1976 == 0 -9.05675 13.81029 -0.656 0.514504
## 1981 - 1976 == 0 -0.78111 18.35985 -0.043 0.966209
## 1982 - 1976 == 0 -2.38339 12.41250 -0.192 0.848390
## 1983 - 1976 == 0 -38.23832 13.92671 -2.746 0.007992 **
## 1984 - 1976 == 0 -28.95867 12.90860 -2.243 0.028641 *
## 1985 - 1976 == 0 -18.33423 15.08488 -1.215 0.229053
## 1986 - 1976 == 0 -28.40711 12.90860 -2.201 0.031688 *
## 1987 - 1976 == 0 -18.64925 12.08019 -1.544 0.127987
## 1988 - 1976 == 0 -16.89143 13.27936 -1.272 0.208362
## 1989 - 1976 == 0 -21.20734 12.30424 -1.724 0.090021 .
## 1990 - 1976 == 0 -22.36611 18.35985 -1.218 0.227994
## 1991 - 1976 == 0 -7.32560 14.34376 -0.511 0.611455
## 1992 - 1976 == 0 -11.15249 12.45224 -0.896 0.374095
## 1993 - 1976 == 0 -17.78203 12.71582 -1.398 0.167222
## 1994 - 1976 == 0 -8.67170 12.78086 -0.678 0.500113
## 1995 - 1976 == 0 13.57536 13.30638 1.020 0.311792
## 1996 - 1976 == 0 -3.34922 13.30638 -0.252 0.802147
## 1997 - 1976 == 0 -11.59916 12.56356 -0.923 0.359644
## 1998 - 1976 == 0 14.62307 18.62181 0.785 0.435440
## 1978 - 1977 == 0 -23.25541 18.23689 -1.275 0.207242
## 1979 - 1977 == 0 10.11959 18.23689 0.555 0.581063
## 1980 - 1977 == 0 0.40936 17.14223 0.024 0.981029
## 1981 - 1977 == 0 8.68500 20.93138 0.415 0.679700
## 1982 - 1977 == 0 7.08272 16.43022 0.431 0.667982
## 1983 - 1977 == 0 -28.77221 17.29679 -1.663 0.101527
## 1984 - 1977 == 0 -19.49257 16.66532 -1.170 0.246846
## 1985 - 1977 == 0 -8.86813 18.12710 -0.489 0.626500
## 1986 - 1977 == 0 -18.94101 16.66532 -1.137 0.260321
## 1987 - 1977 == 0 -9.18314 15.95407 -0.576 0.567075
## 1988 - 1977 == 0 -7.42533 16.79651 -0.442 0.660051
## 1989 - 1977 == 0 -11.74124 16.16342 -0.726 0.470461
## 1990 - 1977 == 0 -12.90000 20.93138 -0.616 0.540067
## 1991 - 1977 == 0 2.14051 17.69601 0.121 0.904134
## 1992 - 1977 == 0 -1.68638 16.10151 -0.105 0.916942
## 1993 - 1977 == 0 -8.31592 16.53313 -0.503 0.616847
## 1994 - 1977 == 0 0.79441 16.37629 0.049 0.961474
## 1995 - 1977 == 0 23.04147 16.81692 1.370 0.175836
## 1996 - 1977 == 0 6.11688 16.81692 0.364 0.717357
## 1997 - 1977 == 0 -2.13305 16.22072 -0.132 0.895826
## 1998 - 1977 == 0 24.08918 21.30941 1.130 0.262862
## 1979 - 1978 == 0 33.37500 14.80072 2.255 0.027861 *
## 1980 - 1978 == 0 23.66477 13.52755 1.749 0.085426 .
## 1981 - 1978 == 0 31.94041 18.23689 1.751 0.085071 .
## 1982 - 1978 == 0 30.33813 12.59792 2.408 0.019176 *
## 1983 - 1978 == 0 -5.51680 13.52755 -0.408 0.684882
## 1984 - 1978 == 0 3.76284 12.90860 0.291 0.771692
## 1985 - 1978 == 0 14.38729 14.93497 0.963 0.339314
## 1986 - 1978 == 0 4.31440 12.90860 0.334 0.739393
## 1987 - 1978 == 0 14.07227 11.92988 1.180 0.242901
## 1988 - 1978 == 0 15.83008 12.90266 1.227 0.224739
## 1989 - 1978 == 0 11.51417 12.16281 0.947 0.347668
## 1990 - 1978 == 0 10.35541 18.23689 0.568 0.572307
## 1991 - 1978 == 0 25.39592 14.02252 1.811 0.075220 .
## 1992 - 1978 == 0 21.56903 12.12366 1.779 0.080377 .
## 1993 - 1978 == 0 14.93949 12.64762 1.181 0.242259
## 1994 - 1978 == 0 24.04982 12.43222 1.934 0.057855 .
## 1995 - 1978 == 0 46.29688 12.85667 3.601 0.000651 ***
## 1996 - 1978 == 0 29.37229 12.85667 2.285 0.025949 *
## 1997 - 1978 == 0 21.12236 12.22384 1.728 0.089226 .
## 1998 - 1978 == 0 47.34459 18.23689 2.596 0.011882 *
## 1980 - 1979 == 0 -9.71023 13.52755 -0.718 0.475706
## 1981 - 1979 == 0 -1.43459 18.23689 -0.079 0.937566
## 1982 - 1979 == 0 -3.03687 12.59792 -0.241 0.810344
## 1983 - 1979 == 0 -38.89180 13.52755 -2.875 0.005611 **
## 1984 - 1979 == 0 -29.61216 12.90860 -2.294 0.025369 *
## 1985 - 1979 == 0 -18.98771 14.93497 -1.271 0.208590
## 1986 - 1979 == 0 -29.06060 12.90860 -2.251 0.028108 *
## 1987 - 1979 == 0 -19.30273 11.92988 -1.618 0.110993
## 1988 - 1979 == 0 -17.54492 12.90266 -1.360 0.179070
## 1989 - 1979 == 0 -21.86083 12.16281 -1.797 0.077398 .
## 1990 - 1979 == 0 -23.01959 18.23689 -1.262 0.211824
## 1991 - 1979 == 0 -7.97908 14.02252 -0.569 0.571503
## 1992 - 1979 == 0 -11.80597 12.12366 -0.974 0.334134
## 1993 - 1979 == 0 -18.43551 12.64762 -1.458 0.150246
## 1994 - 1979 == 0 -9.32518 12.43222 -0.750 0.456186
## 1995 - 1979 == 0 12.92188 12.85667 1.005 0.318967
## 1996 - 1979 == 0 -4.00271 12.85667 -0.311 0.756645
## 1997 - 1979 == 0 -12.25264 12.22384 -1.002 0.320265
## 1998 - 1979 == 0 13.96959 18.23689 0.766 0.446726
## 1981 - 1980 == 0 8.27564 17.14223 0.483 0.631050
## 1982 - 1980 == 0 6.67336 11.04000 0.604 0.547848
## 1983 - 1980 == 0 -29.18157 12.15792 -2.400 0.019560 *
## 1984 - 1980 == 0 -19.90193 11.39114 -1.747 0.085816 .
## 1985 - 1980 == 0 -9.27748 13.57665 -0.683 0.497066
## 1986 - 1980 == 0 -19.35036 11.39114 -1.699 0.094640 .
## 1987 - 1980 == 0 -9.59250 10.28678 -0.933 0.354875
## 1988 - 1980 == 0 -7.83468 11.45073 -0.684 0.496522
## 1989 - 1980 == 0 -12.15059 10.57357 -1.149 0.255132
## 1990 - 1980 == 0 -13.30936 17.14223 -0.776 0.440610
## 1991 - 1980 == 0 1.73115 12.71073 0.136 0.892130
## 1992 - 1980 == 0 -2.09574 10.51193 -0.199 0.842661
## 1993 - 1980 == 0 -8.72528 11.12874 -0.784 0.436158
## 1994 - 1980 == 0 0.38505 10.88700 0.035 0.971906
## 1995 - 1980 == 0 22.63211 11.42620 1.981 0.052291 .
## 1996 - 1980 == 0 5.70753 11.42620 0.500 0.619276
## 1997 - 1980 == 0 -2.54241 10.64946 -0.239 0.812138
## 1998 - 1980 == 0 23.67982 17.29679 1.369 0.176179
## 1982 - 1981 == 0 -1.60228 16.43022 -0.098 0.922644
## 1983 - 1981 == 0 -37.45721 17.29679 -2.166 0.034399 *
## 1984 - 1981 == 0 -28.17757 16.66532 -1.691 0.096154 .
## 1985 - 1981 == 0 -17.55313 18.12710 -0.968 0.336830
## 1986 - 1981 == 0 -27.62601 16.66532 -1.658 0.102687
## 1987 - 1981 == 0 -17.86814 15.95407 -1.120 0.267263
## 1988 - 1981 == 0 -16.11033 16.79651 -0.959 0.341399
## 1989 - 1981 == 0 -20.42624 16.16342 -1.264 0.211297
## 1990 - 1981 == 0 -21.58500 20.93138 -1.031 0.306643
## 1991 - 1981 == 0 -6.54449 17.69601 -0.370 0.712834
## 1992 - 1981 == 0 -10.37138 16.10151 -0.644 0.521991
## 1993 - 1981 == 0 -17.00092 16.53313 -1.028 0.308009
## 1994 - 1981 == 0 -7.89059 16.37629 -0.482 0.631709
## 1995 - 1981 == 0 14.35647 16.81692 0.854 0.396729
## 1996 - 1981 == 0 -2.56812 16.81692 -0.153 0.879148
## 1997 - 1981 == 0 -10.81805 16.22072 -0.667 0.507418
## 1998 - 1981 == 0 15.40418 21.30941 0.723 0.472609
## 1983 - 1982 == 0 -35.85493 11.10858 -3.228 0.002039 **
## 1984 - 1982 == 0 -26.57529 9.96432 -2.667 0.009861 **
## 1985 - 1982 == 0 -15.95085 12.66577 -1.259 0.212857
## 1986 - 1982 == 0 -26.02373 9.96432 -2.612 0.011408 *
## 1987 - 1982 == 0 -16.26586 8.81152 -1.846 0.069914 .
## 1988 - 1982 == 0 -14.50805 10.30126 -1.408 0.164268
## 1989 - 1982 == 0 -18.82396 9.10569 -2.067 0.043104 *
## 1990 - 1982 == 0 -19.98272 16.43022 -1.216 0.228746
## 1991 - 1982 == 0 -4.94221 11.64016 -0.425 0.672685
## 1992 - 1982 == 0 -8.76910 9.26628 -0.946 0.347832
## 1993 - 1982 == 0 -15.39864 9.67190 -1.592 0.116706
## 1994 - 1982 == 0 -6.28831 9.68212 -0.649 0.518551
## 1995 - 1982 == 0 15.95875 10.29949 1.549 0.126617
## 1996 - 1982 == 0 -0.96584 10.29949 -0.094 0.925606
## 1997 - 1982 == 0 -9.21577 9.39664 -0.981 0.330720
## 1998 - 1982 == 0 17.00646 16.56830 1.026 0.308872
## 1984 - 1983 == 0 9.27964 11.46179 0.810 0.421414
## 1985 - 1983 == 0 19.90409 13.77128 1.445 0.153655
## 1986 - 1983 == 0 9.83121 11.46179 0.858 0.394508
## 1987 - 1983 == 0 19.58907 10.32922 1.896 0.062797 .
## 1988 - 1983 == 0 21.34689 11.38881 1.874 0.065832 .
## 1989 - 1983 == 0 17.03098 10.57989 1.610 0.112791
## 1990 - 1983 == 0 15.87221 17.29679 0.918 0.362544
## 1991 - 1983 == 0 30.91272 12.63195 2.447 0.017394 *
## 1992 - 1983 == 0 27.08583 10.55141 2.567 0.012813 *
## 1993 - 1983 == 0 20.45629 11.13294 1.837 0.071180 .
## 1994 - 1983 == 0 29.56662 10.88394 2.717 0.008643 **
## 1995 - 1983 == 0 51.81368 11.30914 4.582 2.44e-05 ***
## 1996 - 1983 == 0 34.88910 11.30914 3.085 0.003096 **
## 1997 - 1983 == 0 26.63916 10.64424 2.503 0.015116 *
## 1998 - 1983 == 0 52.86139 17.14223 3.084 0.003108 **
## 1985 - 1984 == 0 10.62444 12.96928 0.819 0.415969
## 1986 - 1984 == 0 0.55156 10.46569 0.053 0.958147
## 1987 - 1984 == 0 10.30942 9.32873 1.105 0.273591
## 1988 - 1984 == 0 12.06724 10.69153 1.129 0.263604
## 1989 - 1984 == 0 7.75133 9.62096 0.806 0.423668
## 1990 - 1984 == 0 6.59257 16.66532 0.396 0.693837
## 1991 - 1984 == 0 21.63308 12.00358 1.802 0.076620 .
## 1992 - 1984 == 0 17.80619 9.69120 1.837 0.071195 .
## 1993 - 1984 == 0 11.17665 10.18459 1.097 0.276922
## 1994 - 1984 == 0 20.28698 10.09227 2.010 0.048994 *
## 1995 - 1984 == 0 42.53404 10.68062 3.982 0.000190 ***
## 1996 - 1984 == 0 25.60945 10.68062 2.398 0.019679 *
## 1997 - 1984 == 0 17.35952 9.82504 1.767 0.082423 .
## 1998 - 1984 == 0 43.58175 16.81001 2.593 0.011990 *
## 1986 - 1985 == 0 -10.07288 12.96928 -0.777 0.440455
## 1987 - 1985 == 0 -0.31502 12.04166 -0.026 0.979217
## 1988 - 1985 == 0 1.44280 13.13743 0.110 0.912922
## 1989 - 1985 == 0 -2.87311 12.31769 -0.233 0.816374
## 1990 - 1985 == 0 -4.03187 18.12710 -0.222 0.824753
## 1991 - 1985 == 0 11.00863 14.26948 0.771 0.443501
## 1992 - 1985 == 0 7.18174 12.23634 0.587 0.559497
## 1993 - 1985 == 0 0.55221 12.79898 0.043 0.965732
## 1994 - 1985 == 0 9.66253 12.59573 0.767 0.446065
## 1995 - 1985 == 0 31.90959 13.16352 2.424 0.018430 *
## 1996 - 1985 == 0 14.98501 13.16352 1.138 0.259566
## 1997 - 1985 == 0 6.73508 12.39278 0.543 0.588856
## 1998 - 1985 == 0 32.95730 18.56234 1.775 0.080974 .
## 1987 - 1986 == 0 9.75786 9.32873 1.046 0.299827
## 1988 - 1986 == 0 11.51568 10.69153 1.077 0.285827
## 1989 - 1986 == 0 7.19977 9.62096 0.748 0.457227
## 1990 - 1986 == 0 6.04101 16.66532 0.362 0.718281
## 1991 - 1986 == 0 21.08151 12.00358 1.756 0.084232 .
## 1992 - 1986 == 0 17.25462 9.69120 1.780 0.080153 .
## 1993 - 1986 == 0 10.62509 10.18459 1.043 0.301088
## 1994 - 1986 == 0 19.73541 10.09227 1.955 0.055266 .
## 1995 - 1986 == 0 41.98247 10.68062 3.931 0.000225 ***
## 1996 - 1986 == 0 25.05789 10.68062 2.346 0.022350 *
## 1997 - 1986 == 0 16.80796 9.82504 1.711 0.092385 .
## 1998 - 1986 == 0 43.03019 16.81001 2.560 0.013055 *
## 1988 - 1987 == 0 1.75782 9.36462 0.188 0.851749
## 1989 - 1987 == 0 -2.55809 8.24140 -0.310 0.757354
## 1990 - 1987 == 0 -3.71686 15.95407 -0.233 0.816589
## 1991 - 1987 == 0 11.32365 10.67732 1.061 0.293225
## 1992 - 1987 == 0 7.49676 8.27310 0.906 0.368536
## 1993 - 1987 == 0 0.86722 8.79756 0.099 0.921809
## 1994 - 1987 == 0 9.97755 8.71342 1.145 0.256800
## 1995 - 1987 == 0 32.22461 9.44515 3.412 0.001171 **
## 1996 - 1987 == 0 15.30003 9.44515 1.620 0.110591
## 1997 - 1987 == 0 7.05009 8.33792 0.846 0.401222
## 1998 - 1987 == 0 33.27232 16.03612 2.075 0.042372 *
## 1989 - 1988 == 0 -4.31591 9.62361 -0.448 0.655456
## 1990 - 1988 == 0 -5.47467 16.79651 -0.326 0.745622
## 1991 - 1988 == 0 9.56583 11.53751 0.829 0.410386
## 1992 - 1988 == 0 5.73894 9.60194 0.598 0.552337
## 1993 - 1988 == 0 -0.89059 10.05120 -0.089 0.929695
## 1994 - 1988 == 0 8.21973 9.94532 0.826 0.411855
## 1995 - 1988 == 0 30.46679 10.56650 2.883 0.005483 **
## 1996 - 1988 == 0 13.54221 10.56650 1.282 0.204991
## 1997 - 1988 == 0 5.29228 9.57737 0.553 0.582637
## 1998 - 1988 == 0 31.51450 16.66974 1.891 0.063602 .
## 1990 - 1989 == 0 -1.15876 16.16342 -0.072 0.943091
## 1991 - 1989 == 0 13.88175 10.86908 1.277 0.206543
## 1992 - 1989 == 0 10.05485 8.60023 1.169 0.247050
## 1993 - 1989 == 0 3.42532 9.05878 0.378 0.706697
## 1994 - 1989 == 0 12.53564 9.01126 1.391 0.169416
## 1995 - 1989 == 0 34.78271 9.70911 3.582 0.000690 ***
## 1996 - 1989 == 0 17.85812 9.70911 1.839 0.070902 .
## 1997 - 1989 == 0 9.60819 8.63586 1.113 0.270396
## 1998 - 1989 == 0 35.83042 16.17582 2.215 0.030629 *
## 1991 - 1990 == 0 15.04051 17.69601 0.850 0.398796
## 1992 - 1990 == 0 11.21362 16.10151 0.696 0.488893
## 1993 - 1990 == 0 4.58408 16.53313 0.277 0.782545
## 1994 - 1990 == 0 13.69441 16.37629 0.836 0.406398
## 1995 - 1990 == 0 35.94147 16.81692 2.137 0.036736 *
## 1996 - 1990 == 0 19.01688 16.81692 1.131 0.262707
## 1997 - 1990 == 0 10.76695 16.22072 0.664 0.509418
## 1998 - 1990 == 0 36.98918 21.30941 1.736 0.087817 .
## 1992 - 1991 == 0 -3.82689 10.87482 -0.352 0.726164
## 1993 - 1991 == 0 -10.45643 10.97923 -0.952 0.344788
## 1994 - 1991 == 0 -1.34610 11.13470 -0.121 0.904187
## 1995 - 1991 == 0 20.90096 11.88976 1.758 0.083952 .
## 1996 - 1991 == 0 3.97638 11.88976 0.334 0.739236
## 1997 - 1991 == 0 -4.27356 10.64946 -0.401 0.689654
## 1998 - 1991 == 0 21.94867 17.52595 1.252 0.215383
## 1993 - 1992 == 0 -6.62954 9.14341 -0.725 0.471280
## 1994 - 1992 == 0 2.48079 8.97050 0.277 0.783092
## 1995 - 1992 == 0 24.72785 9.68708 2.553 0.013298 *
## 1996 - 1992 == 0 7.80327 9.68708 0.806 0.423747
## 1997 - 1992 == 0 -0.44667 8.60178 -0.052 0.958762
## 1998 - 1992 == 0 25.77556 16.17879 1.593 0.116466
## 1994 - 1993 == 0 9.11033 9.50529 0.958 0.341748
## 1995 - 1993 == 0 31.35739 10.30850 3.042 0.003505 **
## 1996 - 1993 == 0 14.43280 10.30850 1.400 0.166724
## 1997 - 1993 == 0 6.18287 9.04464 0.684 0.496907
## 1998 - 1993 == 0 32.40510 16.54163 1.959 0.054844 .
## 1995 - 1994 == 0 22.24706 10.03707 2.216 0.030526 *
## 1996 - 1994 == 0 5.32248 10.03707 0.530 0.597907
## 1997 - 1994 == 0 -2.92746 8.99241 -0.326 0.745919
## 1998 - 1994 == 0 23.29477 16.37018 1.423 0.160002
## 1996 - 1995 == 0 -16.92458 10.46569 -1.617 0.111181
## 1997 - 1995 == 0 -25.17452 9.77607 -2.575 0.012548 *
## 1998 - 1995 == 0 1.04771 16.57783 0.063 0.949821
## 1997 - 1996 == 0 -8.24993 9.77607 -0.844 0.402139
## 1998 - 1996 == 0 17.97229 16.57783 1.084 0.282724
## 1998 - 1997 == 0 26.22223 16.21043 1.618 0.111080
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Univariate p values reported)
ENSO.83_97$Year <- as.factor(ENSO.83_97$Year)
comp_years_ENSO.83_97 <- aov(Coral_cover ~ Year + Thermal_regime, data = ENSO.83_97)
m.comp_ENSO.83_97 <- glht(comp_years_ENSO.83_97, linfct = mcp(Year = "Tukey"))
summary(m.comp_ENSO.83_97, test = univariate())
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = Coral_cover ~ Year + Thermal_regime, data = ENSO.83_97)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 1984 - 1983 == 0 12.351578 10.375355 1.190 0.239968
## 1985 - 1983 == 0 21.853429 12.451269 1.755 0.085898 .
## 1986 - 1983 == 0 12.903140 10.375355 1.244 0.219935
## 1987 - 1983 == 0 21.659577 9.319358 2.324 0.024588 *
## 1988 - 1983 == 0 21.654675 10.234189 2.116 0.039796 *
## 1989 - 1983 == 0 19.121678 9.546497 2.003 0.051090 .
## 1990 - 1983 == 0 17.821554 15.603472 1.142 0.259301
## 1991 - 1983 == 0 31.408598 11.354258 2.766 0.008138 **
## 1992 - 1983 == 0 27.940801 9.494612 2.943 0.005082 **
## 1993 - 1983 == 0 23.016432 10.063858 2.287 0.026843 *
## 1994 - 1983 == 0 30.202718 9.785636 3.086 0.003426 **
## 1995 - 1983 == 0 51.570014 10.163303 5.074 6.85e-06 ***
## 1996 - 1983 == 0 34.645431 10.163303 3.409 0.001366 **
## 1997 - 1983 == 0 27.374437 9.570674 2.860 0.006346 **
## 1998 - 1983 == 0 51.886723 15.418676 3.365 0.001551 **
## 1985 - 1984 == 0 9.501851 11.732643 0.810 0.422188
## 1986 - 1984 == 0 0.551563 9.403948 0.059 0.953483
## 1987 - 1984 == 0 9.307999 8.393339 1.109 0.273206
## 1988 - 1984 == 0 9.303098 9.673243 0.962 0.341211
## 1989 - 1984 == 0 6.770100 8.654255 0.782 0.438055
## 1990 - 1984 == 0 5.469976 15.036262 0.364 0.717685
## 1991 - 1984 == 0 19.057020 10.836676 1.759 0.085303 .
## 1992 - 1984 == 0 15.589224 8.769831 1.778 0.082081 .
## 1993 - 1984 == 0 10.664854 9.156204 1.165 0.250118
## 1994 - 1984 == 0 17.851140 9.129581 1.955 0.056638 .
## 1995 - 1984 == 0 39.218436 9.691117 4.047 0.000197 ***
## 1996 - 1984 == 0 22.293853 9.691117 2.300 0.026008 *
## 1997 - 1984 == 0 15.022859 8.885050 1.691 0.097641 .
## 1998 - 1984 == 0 39.535145 15.199628 2.601 0.012458 *
## 1986 - 1985 == 0 -8.950288 11.732643 -0.763 0.449445
## 1987 - 1985 == 0 -0.193851 10.880330 -0.018 0.985862
## 1988 - 1985 == 0 -0.198753 11.878661 -0.017 0.986723
## 1989 - 1985 == 0 -2.731750 11.148223 -0.245 0.807516
## 1990 - 1985 == 0 -4.031875 16.288116 -0.248 0.805596
## 1991 - 1985 == 0 9.555169 12.906122 0.740 0.462844
## 1992 - 1985 == 0 6.087373 11.030341 0.552 0.583705
## 1993 - 1985 == 0 1.163004 11.599831 0.100 0.920573
## 1994 - 1985 == 0 8.349289 11.367353 0.734 0.466374
## 1995 - 1985 == 0 29.716586 11.930022 2.491 0.016413 *
## 1996 - 1985 == 0 12.792002 11.930022 1.072 0.289200
## 1997 - 1985 == 0 5.521009 11.187549 0.493 0.624010
## 1998 - 1985 == 0 30.033295 16.807767 1.787 0.080549 .
## 1987 - 1986 == 0 8.756437 8.393339 1.043 0.302280
## 1988 - 1986 == 0 8.751535 9.673243 0.905 0.370330
## 1989 - 1986 == 0 6.218538 8.654255 0.719 0.476051
## 1990 - 1986 == 0 4.918413 15.036262 0.327 0.745074
## 1991 - 1986 == 0 18.505457 10.836676 1.708 0.094439 .
## 1992 - 1986 == 0 15.037661 8.769831 1.715 0.093130 .
## 1993 - 1986 == 0 10.113292 9.156204 1.105 0.275108
## 1994 - 1986 == 0 17.299577 9.129581 1.895 0.064402 .
## 1995 - 1986 == 0 38.666874 9.691117 3.990 0.000235 ***
## 1996 - 1986 == 0 21.742290 9.691117 2.244 0.029720 *
## 1997 - 1986 == 0 14.471297 8.885050 1.629 0.110202
## 1998 - 1986 == 0 38.983583 15.199628 2.565 0.013653 *
## 1988 - 1987 == 0 -0.004902 8.444989 -0.001 0.999539
## 1989 - 1987 == 0 -2.537899 7.407768 -0.343 0.733460
## 1990 - 1987 == 0 -3.838024 14.381093 -0.267 0.790755
## 1991 - 1987 == 0 9.749021 9.617622 1.014 0.316046
## 1992 - 1987 == 0 6.281224 7.459867 0.842 0.404143
## 1993 - 1987 == 0 1.356855 7.915355 0.171 0.864645
## 1994 - 1987 == 0 8.543141 7.854396 1.088 0.282399
## 1995 - 1987 == 0 29.910437 8.539886 3.502 0.001038 **
## 1996 - 1987 == 0 12.985854 8.539886 1.521 0.135202
## 1997 - 1987 == 0 5.714860 7.513916 0.761 0.450795
## 1998 - 1987 == 0 30.227146 14.473044 2.089 0.042314 *
## 1989 - 1988 == 0 -2.532997 8.679641 -0.292 0.771727
## 1990 - 1988 == 0 -3.833122 15.150473 -0.253 0.801393
## 1991 - 1988 == 0 9.753922 10.369343 0.941 0.351798
## 1992 - 1988 == 0 6.286126 8.639110 0.728 0.470524
## 1993 - 1988 == 0 1.361757 9.084485 0.150 0.881499
## 1994 - 1988 == 0 8.548042 8.940295 0.956 0.344007
## 1995 - 1988 == 0 29.915339 9.498092 3.150 0.002871 **
## 1996 - 1988 == 0 12.990755 9.498092 1.368 0.178044
## 1997 - 1988 == 0 5.719762 8.609504 0.664 0.509780
## 1998 - 1988 == 0 30.232048 14.998026 2.016 0.049691 *
## 1990 - 1989 == 0 -1.300125 14.584825 -0.089 0.929356
## 1991 - 1989 == 0 12.286919 9.787894 1.255 0.215703
## 1992 - 1989 == 0 8.819123 7.765025 1.136 0.261944
## 1993 - 1989 == 0 3.894754 8.143592 0.478 0.634730
## 1994 - 1989 == 0 11.081039 8.129358 1.363 0.179490
## 1995 - 1989 == 0 32.448336 8.775695 3.698 0.000579 ***
## 1996 - 1989 == 0 15.523753 8.775695 1.769 0.083532 .
## 1997 - 1989 == 0 8.252759 7.788304 1.060 0.294843
## 1998 - 1989 == 0 32.765045 14.592673 2.245 0.029597 *
## 1991 - 1990 == 0 13.587044 15.968790 0.851 0.399260
## 1992 - 1990 == 0 10.119248 14.494919 0.698 0.488614
## 1993 - 1990 == 0 5.194879 14.932861 0.348 0.729515
## 1994 - 1990 == 0 12.381164 14.752999 0.839 0.405680
## 1995 - 1990 == 0 33.748461 15.190776 2.222 0.031268 *
## 1996 - 1990 == 0 16.823877 15.190776 1.108 0.273833
## 1997 - 1990 == 0 9.552884 14.614907 0.654 0.516598
## 1998 - 1990 == 0 34.065170 19.259680 1.769 0.083569 .
## 1992 - 1991 == 0 -3.467796 9.790236 -0.354 0.724800
## 1993 - 1991 == 0 -8.392165 9.902026 -0.848 0.401094
## 1994 - 1991 == 0 -1.205880 10.014648 -0.120 0.904681
## 1995 - 1991 == 0 20.161416 10.688016 1.886 0.065569 .
## 1996 - 1991 == 0 3.236833 10.688016 0.303 0.763372
## 1997 - 1991 == 0 -4.034160 9.577556 -0.421 0.675564
## 1998 - 1991 == 0 20.478126 15.763214 1.299 0.200381
## 1993 - 1992 == 0 -4.924369 8.277348 -0.595 0.554812
## 1994 - 1992 == 0 2.261916 8.062383 0.281 0.780312
## 1995 - 1992 == 0 23.629213 8.730875 2.706 0.009511 **
## 1996 - 1992 == 0 6.704629 8.730875 0.768 0.446458
## 1997 - 1992 == 0 -0.566364 7.731752 -0.073 0.941923
## 1998 - 1992 == 0 23.945922 14.586686 1.642 0.107487
## 1994 - 1993 == 0 7.186285 8.596189 0.836 0.407485
## 1995 - 1993 == 0 28.553582 9.335014 3.059 0.003699 **
## 1996 - 1993 == 0 11.628999 9.335014 1.246 0.219168
## 1997 - 1993 == 0 4.358005 8.178387 0.533 0.596690
## 1998 - 1993 == 0 28.870291 14.934120 1.933 0.059383 .
## 1995 - 1994 == 0 21.367296 9.033108 2.365 0.022280 *
## 1996 - 1994 == 0 4.442713 9.033108 0.492 0.625182
## 1997 - 1994 == 0 -2.828280 8.080350 -0.350 0.727922
## 1998 - 1994 == 0 21.684006 14.744958 1.471 0.148206
## 1996 - 1995 == 0 -16.924583 9.403948 -1.800 0.078464 .
## 1997 - 1995 == 0 -24.195577 8.799034 -2.750 0.008496 **
## 1998 - 1995 == 0 0.316709 14.905050 0.021 0.983139
## 1997 - 1996 == 0 -7.270994 8.799034 -0.826 0.412877
## 1998 - 1996 == 0 17.241292 14.905050 1.157 0.253350
## 1998 - 1997 == 0 24.512286 14.601377 1.679 0.099979 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Univariate p values reported)
Figure1<-grid.arrange(ENSO, Model1_plot, arrangeGrob(
Model2_plot, Model3b_plot, Model4b_plot, widths = c(3.9/12, 3.80/12, 4.30/12), ncol = 3),
nrow=3)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
#ggsave(file="Fig1_V5.svg", plot=Figure1, width=7, height=7)
# Load packages
library(tidyverse)
library(dplyr)
library(data.table)
library(RColorBrewer)
# Select DHW data file: maxDHW.csv
dhw <- read.csv("Data/maxDHW.csv", header=TRUE)
# Exclude years 1981, and 2017-2019
dhw <- dhw[!dhw$Year == 1981, ]
# Create groups for plot
ts.Cano_island <- subset(dhw, dhw$Site == "Cano_island")
ts.Uva_island <- subset(dhw, dhw$Site == "Uva_Island")
ts.Gorgona_island <- subset(dhw, dhw$Site == "Gorgona_Island")
ts.Darwin_island <- subset(dhw, dhw$Site == "Darwin_North_Anchorage")
ts.Bahia_Banderas <- subset(dhw, dhw$Site == "Bahia_Banderas")
# Year, maxDHW
par(oma = c(1,1,1,1))
par(mar = c(4,4,1,1), cex.axis = 1, las = 1)
par(lwd = 1)
plot(ts.Cano_island$Year,ts.Cano_island$maxDHW, # Create a white false line with real data
xlab="Years", ylab="Degree heating weeks (°C-weeks)", cex=1.5, cex.lab=0.6, col="white",axes = FALSE,
ylim = c(0, 28), lwd = 1.5)
par(new=T)
rect(1982, 0, 1983.5, 28, col="#fee8c8", border ="#fee8c8") # Strong ENSO years
par(new=T)
rect(1997, 0, 1998.5, 28, col="#fee8c8", border ="#fee8c8")
par(new=T)
plot(ts.Cano_island$Year, ts.Cano_island$maxDHW,type="l", # Thermally stable
xlab="", ylab="", lwd = 1.5, col="#8c2d04",axes = FALSE,
ylim = c(0, 28))
par(new=T)
plot(ts.Uva_island$Year, ts.Uva_island$maxDHW, type="l", # Thermally stable
xlab="", ylab="", lwd = 1.5, col="#feb24c",axes = FALSE,
ylim = c(0, 28))
par(new=T)
plot(ts.Gorgona_island$Year,ts.Gorgona_island$maxDHW, type="l", # Tropical upwelling
xlab="", ylab="",lwd = 1.5, col="#fb6a4a", axes = FALSE,
ylim = c(0, 28))
par(new=T)
plot(ts.Darwin_island$Year,ts.Darwin_island$maxDHW, type="l", # Galapagos
xlab="", ylab="", lwd = 1.5, col="#9e9ac8",axes = FALSE,
ylim = c(0, 28))
par(new=T)
plot(ts.Bahia_Banderas$Year, ts.Bahia_Banderas$maxDHW, type="l", # Seasonal
xlab="", ylab="",lwd = 1.5,col="#a6d96a", axes = FALSE,
ylim = c(0, 28))
y <- c(0,4,8,12,16,20,24,28) # Axis
axis(2,lwd = 0.5, at = y,cex.axis=0.6)
box(lwd = 0.5)
x <-c(1982:2014)
YearsR = as.character(c(1982:2014))
axis(1, at = x, labels = YearsR, lwd = 0.4, cex.axis=0.5, las = 2)
abline(h=4, col="#525252", lty=2)
abline(h=8, col="#525252", lty=2)
# Legend
sites <- c("Thermally stable - Caño Island","Thermally stable - Uva Island",
"Tropical upwelling - Gorgona Island","Equatorial Upwelling - Darwin Island",
"Seasonal - BahÃa Banderas")
par(xpd=TRUE)
legend(x = 2000, y = 20,
xpd = NA, # or TRUE
inset=c(-0.5,0),
cex = 0.4,
legend = sites, # box categories
col = c("#8c2d04","#feb24c","#fb6a4a","#9e9ac8","#a6d96a"),
lty= 1,
lwd = 1.5,
x.intersp = 0.5,
y.intersp = 2,
bty="n")
dev.off() # reset parameters for next plot
## null device
## 1
We tested the hypothesis that the strength of the Log_annual_rate_of_change (response variable) is a function of the maxDHW experienced. We had measured the rate of change from four Thermal regimes, and fitted Thermal regimes as a random intercept, and estimate a common slope (change in the rate of change) for maxDHW across all Thermal regimes by fitting it as a fixed effect.
## Load data
roc.cover <- read.csv("Data/RoC_DHW.csv", header=TRUE) ### Select file: RoC_DWH.csv
# Load packages
library(lme4)
library(lmerTest)
library(MuMIn)
model6 <- lmerTest::lmer(Log_annual_rate_of_change ~ maxDHW + (1|Thermal_regime), data=roc.cover)
step(model6)
## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 4 -17.738 43.476
## (1 | Thermal_regime) 0 3 -27.772 61.545 20.069 1 7.469e-06
##
## <none>
## (1 | Thermal_regime) ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## maxDHW 0 0.65054 0.65054 1 128.59 9.873 0.002082 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## Log_annual_rate_of_change ~ maxDHW + (1 | Thermal_regime)
summary(model6)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Log_annual_rate_of_change ~ maxDHW + (1 | Thermal_regime)
## Data: roc.cover
##
## REML criterion at convergence: 35.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.7962 -0.1472 0.1167 0.4152 1.8621
##
## Random effects:
## Groups Name Variance Std.Dev.
## Thermal_regime (Intercept) 0.02501 0.1581
## Residual 0.06589 0.2567
## Number of obs: 133, groups: Thermal_regime, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.085763 0.084256 3.155381 -1.018 0.38032
## maxDHW -0.016316 0.005193 128.594887 -3.142 0.00208 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## maxDHW -0.165
r.squaredGLMM(model6) # returns the marginal and the conditional R2
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help
## page.
## R2m R2c
## [1,] 0.05249512 0.3131505
# marginal R2 describes the proportion of variance explained by the fixed factor(s) alone:
# conditional R2 describes the proportion of variance explained by both the fixed and random factors
# Nakagawa, S. & Schielzeth, H. (2013).
# Load packages
library(tidyverse)
library(dplyr)
library(data.table)
library(RColorBrewer)
# Create groups for plot
Seasonal = subset(roc.cover, roc.cover$Thermal_regime == "Seasonal")
Thermally_Stable = subset(roc.cover, roc.cover$Thermal_regime == "Thermally_Stable")
Tropical_Upwelling = subset(roc.cover, roc.cover$Thermal_regime == "Upwelling")
Galapagos = subset(roc.cover, roc.cover$Thermal_regime == "Galapagos_Islands")
# Plot
par(oma = c(1,1,1,1))
par(mar = c(4,4,1,1), cex.axis = 1, las = 1)
par(lwd = 1)
plot (Seasonal$maxDHW, Seasonal$Log_annual_rate_of_change,
pch = 21, lwd = 0.5, bg = "#a6d96a93", col="black", axes = FALSE,
cex=1.5,cex.lab=0.6,
xlim = c(0,30), ylim = c(-1.75, 1),
xlab = "Degree heating weeks (°C-weeks)",
ylab = "Annual rate of change in coral cover")
par(new=T);
plot (Thermally_Stable$maxDHW, Thermally_Stable$Log_annual_rate_of_change,
xlab="", ylab="", pch = 21, lwd = 0.5, col="black", bg = "#feb24c93", axes = FALSE,
cex=1.5,cex.lab=0.6,
xlim = c(0,30), ylim = c(-1.75, 1))
par(new=T);
plot (Tropical_Upwelling$maxDHW, Tropical_Upwelling$Log_annual_rate_of_change,
xlab="", ylab="", pch = 21, lwd = 0.5, bg = "#fb6a4a93",col="black", axes = FALSE,
cex=1.5,cex.lab=0.6,
xlim = c(0,30), ylim = c(-1.75, 1))
par(new=T);
plot (Galapagos$maxDHW, Galapagos$Log_annual_rate_of_change,
xlab="", ylab="", pch = 21, lwd = 0.5, bg = "#9e9ac893",col="black", axes = FALSE,
cex=1.5,cex.lab=0.6,
xlim = c(0,30), ylim = c(-1.75, 1))
box(lwd = 0.5)
x <-c(0,5,10,15,20,25,30)
axis(1, lwd = 0.5, at = x, cex.axis = 0.5)
y <- c(-1.5,-1.0,-0.5,0,0.5,1)
axis(2,lwd = 0.5, at = y, cex.axis = 0.5)
abline(v=4, col="#525252", lty=2); abline(v=8, col="#525252", lty=2); abline(h=0, col="#525252", lty=2);
# Axis text I, II, III, IV
text(3.5, 0.3, "II", cex = 0.6); text(5, 0.3, "I", cex = 0.6)
text(3.5, - 0.3, "III", cex = 0.6); text(5, -0.3, "IV", cex = 0.6)
par(xpd=TRUE)
legend(x = 10, y = -0.5,
xpd = NA,
inset=c(-0.5,0),
legend = c("Equatorial Upwelling", "Thermally stable", "Tropical upwelling","Seasonal"),
pch=21,
pt.bg = c("#9e9ac8","#fb6a4a","#fb6a4a","#a6d96a"),
x.intersp = 0.5,
y.intersp = 2,
bty="n")
dev.off() # reset parameters for next plot
## null device
## 1
No stress
no.stress <- subset(roc.cover, maxDHW < 4) #subset of DHW < 4
length(no.stress$maxDHW)/length(roc.cover$maxDHW) # ratio no stress: total
## [1] 0.8496241
1- length(no.stress$maxDHW)/length(roc.cover$maxDHW) # ratio stress: total
## [1] 0.1503759
positive.rate <- subset(no.stress, Log_annual_rate_of_change > 0) #subset of Rate > 0
length(positive.rate$Log_annual_rate_of_change) # Number of observations in
## [1] 48
negative.rate <- subset(no.stress, Log_annual_rate_of_change < 0) #subset of Rate > 0
length(negative.rate$Log_annual_rate_of_change) # Number of observations in
## [1] 58
No stress seasonal
no.stress.seasonal <- subset(no.stress, Thermal_regime == "Seasonal")
positive.rate <- subset(no.stress.seasonal, Log_annual_rate_of_change > 0)
negative.rate <- subset(no.stress.seasonal, Log_annual_rate_of_change < 0)
length(positive.rate$Log_annual_rate_of_change)/length(negative.rate$Log_annual_rate_of_change)
## [1] 0.07142857
No stress upwelling
no.stress.upwelling <- subset(no.stress, Thermal_regime == "Upwelling")
positive.rate <- subset(no.stress.upwelling, Log_annual_rate_of_change > 0)
negative.rate <- subset(no.stress.upwelling, Log_annual_rate_of_change < 0)
length(positive.rate$Log_annual_rate_of_change)/length(negative.rate$Log_annual_rate_of_change)
## [1] 1
No stress thermally stable
no.stress.Thermally_Stable <- subset(no.stress, Thermal_regime == "Thermally_Stable")
positive.rate <- subset(no.stress.Thermally_Stable, Log_annual_rate_of_change > 0)
negative.rate <- subset(no.stress.Thermally_Stable, Log_annual_rate_of_change < 0)
length(positive.rate$Log_annual_rate_of_change)/length(negative.rate$Log_annual_rate_of_change)
## [1] 1.2
No stress Galapagos
no.stress.Galapagos_Islands <- subset(no.stress, Thermal_regime == "Galapagos_Islands")
positive.rate <- subset(no.stress.Galapagos_Islands, Log_annual_rate_of_change > 0)
negative.rate <- subset(no.stress.Galapagos_Islands, Log_annual_rate_of_change < 0)
length(positive.rate$Log_annual_rate_of_change)/length(negative.rate$Log_annual_rate_of_change)
## [1] 1
library(pals)
library(dplyr)
library(reshape2)
library(RColorBrewer)
# Barplot number of surveys per year
obs_num = as.data.frame(table(cover$Year,cover$Region))
colnames(obs_num) = c("Year", "Region", "Freq")
wide_obs_num = dcast(obs_num, Region~Year, value = "Freq")
## Using Freq as value column: use value.var to override.
df = data.matrix(wide_obs_num[,2:45])
as.table(df)
## 1970 1971 1972 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984
## A 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## B 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## C 0 0 1 0 0 0 0 0 1 1 0 0 1 0
## D 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## E 0 0 0 0 0 0 0 0 0 1 0 1 0 4
## F 0 0 0 2 3 2 0 0 0 0 0 1 0 1
## G 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## H 3 10 2 13 4 4 3 5 2 5 2 3 5 3
## I 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## J 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## K 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## L 0 0 0 0 0 0 0 2 0 0 0 0 0 0
## M 0 0 0 0 0 0 0 0 0 0 0 1 0 0
## N 2 6 0 2 1 0 0 0 0 0 0 1 4 2
## O 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## P 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998
## A 0 0 0 0 0 0 0 0 0 1 0 0 0 0
## B 0 0 5 0 0 0 0 1 0 0 0 0 0 0
## C 0 0 1 2 4 0 0 1 0 0 2 2 0 4
## D 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## E 12 1 2 0 1 0 0 11 0 2 0 6 0 1
## F 0 1 1 0 1 0 0 0 2 0 0 0 0 0
## G 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## H 4 5 2 3 1 1 0 4 1 4 4 0 1 1
## I 0 0 0 0 0 0 2 0 0 0 0 0 16 10
## J 0 0 1 1 1 0 3 1 1 0 0 0 2 2
## K 0 0 0 0 0 0 0 0 0 0 0 0 2 1
## L 0 0 0 0 0 0 0 0 0 1 1 4 4 2
## M 0 0 0 0 0 0 0 1 0 0 0 0 0 0
## N 0 1 5 1 2 0 1 1 1 1 1 1 1 1
## O 0 0 0 0 0 0 0 0 1 1 0 0 0 0
## P 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012
## A 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## B 0 0 0 1 0 0 0 1 0 0 4 0 0 0
## C 7 4 3 4 6 4 3 6 3 2 3 2 2 3
## D 0 0 0 0 0 0 0 0 0 0 0 4 0 1
## E 1 0 2 0 4 2 2 0 1 0 2 0 0 0
## F 0 0 1 1 0 1 3 6 2 2 0 0 1 0
## G 0 0 0 0 0 0 0 0 0 0 1 0 0 0
## H 0 1 4 14 12 1 1 1 1 0 1 2 0 0
## I 0 0 0 6 3 0 1 1 4 0 17 6 2 0
## J 0 0 0 1 2 0 0 0 0 5 7 0 0 1
## K 1 5 5 3 2 0 0 0 0 0 2 0 0 0
## L 0 0 2 0 1 0 1 3 2 0 17 3 1 1
## M 0 0 0 0 0 1 2 14 8 0 0 0 0 1
## N 1 1 1 1 1 1 0 0 0 0 1 0 0 0
## O 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## P 0 0 0 0 0 2 0 0 0 0 0 0 0 0
## 2013 2014
## A 0 0
## B 1 1
## C 2 1
## D 4 2
## E 0 0
## F 0 0
## G 0 0
## H 0 0
## I 2 2
## J 0 0
## K 0 0
## L 1 3
## M 0 0
## N 0 0
## O 0 0
## P 0 0
df[ df>1 ] <- 1
levels(cover$Region) # How many regions? = how many colours
## [1] "Clipperton" "Cocos_Islands"
## [3] "Colombia" "Costa_Rica_Central"
## [5] "Costa_Rica_Southern" "Eastern_Galapagos_Islands"
## [7] "El_Salvador" "Gulf_of_Chiriqui"
## [9] "Mexican_Central" "Mexican_Northern"
## [11] "Mexican_Southern" "Nicaragua_Papagayo_Zone"
## [13] "Northern_Galapagos_Islands" "Panama_Bight"
## [15] "Revillagigedo_Islands" "Western_Galapagos_Islands"
Regions <- c(levels(obs_num$Region))
color.regions = warmcool(17)
# Figure 2a
par(mfcol=c(2,1), cex.main = 1.5, mar = c(4,4,1,1), cex.lab = 1.5, cex.axis = 1, las = 1)
barplot(df,
ylim=c(0,15),
col = color.regions,
ylab="Number of regions surveyed per year",
xlab="Year",
#lwd = 1,
las = 2,
axis.lty = 1,
cex.names = 0.4,
cex.axis = 0.8)
# Legend as a box
legend( #"bottomright",
x = 2, y = 16,
xpd = NA, # or TRUE
legend = Regions,
fill = color.regions[1:17],
inset=c(-0.5,0),
cex = 0.35,
text.font=0.5,
x.intersp = -0.5,
y.intersp = 1,
lwd = 0.1,
bty="n")
# Figure 2b
obs_num = as.data.frame(table(cover$Year,cover$Country)) # Table per country
colnames(obs_num) = c("Year", "Country", "Freq")
wide_obs_num = dcast(obs_num, Country~Year, value = "Freq")
## Using Freq as value column: use value.var to override.
dfc = data.matrix(wide_obs_num[,2:45])
as.table(dfc)
## 1970 1971 1972 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984
## A 0 0 1 0 0 0 0 0 1 1 0 0 1 0
## B 0 0 2 9 0 0 0 2 0 2 0 1 0 4
## C 0 0 0 2 3 2 0 0 0 0 0 2 0 1
## D 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## E 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## F 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## G 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## H 5 16 0 6 5 4 3 5 2 4 2 4 9 5
## 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998
## A 0 0 1 2 4 0 0 1 0 0 2 2 0 4
## B 12 1 7 0 1 0 0 15 0 6 4 10 4 3
## C 0 1 1 0 1 0 0 1 2 0 0 0 0 0
## D 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## E 0 0 0 0 0 0 0 0 0 1 0 0 0 0
## F 0 0 1 1 1 0 5 1 2 1 0 0 20 13
## G 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## H 4 6 7 4 3 1 1 2 2 2 2 1 2 2
## 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012
## A 7 4 3 4 6 4 3 6 3 2 3 2 2 3
## B 1 0 7 1 5 2 3 4 3 0 12 7 1 2
## C 0 0 1 1 0 4 5 20 10 2 0 0 1 1
## D 0 0 0 0 0 0 0 0 0 0 1 0 0 0
## E 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## F 1 5 5 10 7 0 1 1 4 5 26 6 2 1
## G 0 0 0 0 0 0 0 0 0 0 11 0 0 0
## H 1 2 2 15 13 2 1 1 1 0 2 2 0 0
## 2013 2014
## A 2 1
## B 6 6
## C 0 0
## D 0 0
## E 0 0
## F 2 2
## G 0 0
## H 0 0
dfc[ dfc>1 ] <- 1 # replace all values greater than one with one
color.countries <- brewer.pal(9, "Oranges")
countries = c(levels(obs_num$Country))
barplot(dfc,
col = color.countries,
ylim=c(0,10),
ylab="Number of countries surveyed per year",
xlab="Year",
las = 2,
axis.lty = 1,
cex.names = 0.4,
cex.axis = 0.8)
# Legend as a box
legend( #"bottomright",
x = 2, y = 10,
xpd = NA, # or TRUE
legend = countries,
fill = color.countries[1:9],
inset=c(-0.5,0),
cex = 0.35,
text.font=1,
x.intersp = -0.5,
y.intersp = 1,
lwd = 0.1,
bty="n")
dev.off() # reset parameters for next plot
## null device
## 1
Temporal variation in the live coral cover for 1970–2014. * (a) Tropical upwelling regime. * (b) Seasonal regime. * (c) Thermally stable regime. * (d) Equatorial upwelling regime. The smooth line represents the best fit to the time series with a 95% confidence level interval and a span = 0.3.
library(ggplot2)
library(lemon)
aggr.region$Thermal_regime<-factor(aggr.region$Thermal_regime,
levels=c("Tropical upwelling", "Seasonal",
"Thermally stable", "Equatorial upwelling"),
labels = c ("a", "b", "c", "d"))
summary(aggr.region$Thermal_regime)
## a b c d
## 75 28 74 25
FigureS3<-ggplot(aggr.region, aes(x=Year, y=Coral_cover)) +
theme_bw() + theme(panel.border = element_blank(), # set ggplot to white background
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(),
strip.background = element_blank(),
strip.text = element_text(face="bold",hjust = 0, vjust = -1))+
geom_point()+ geom_smooth (span = 0.3) +
scale_x_continuous(limits=c(1970, 2014),
breaks = seq(1970, 2014, by=5),
expand = c(0.01, 0.01), "Year") +
scale_y_continuous(limits=c(-19, 91), expand = c(0, 0), "Live coral cover (%)") +
facet_wrap(Thermal_regime~., ncol=1, scales = "free_x", strip.position ="top")
FigureS3
#ggsave(file="FigureS3.svg", plot=FigureS3, width=6.69, height=6.69)